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Formation of sprays from conical liquid sheets 

By Bill Peck AND N. N. Mansour 
1. Motivation and objectives . 

Our objective is to predict droplet size distributions created by fuel injector noz- 
zles in jet turbines. These results will be used to determine the initial conditions 
for numerical simulations of the combustion process in gas turbine combustors. 
To predict the droplet size distribution, we are currently constructing a numerical 
model to understand the instability and breakup of thin conical liquid sheets. This 
geometry serves as a simplified model of the liquid jet emerging from a real nozzle. 

The physics of this process is difficult to study experimentally as the time and 
length scales are very short. From existing photographic data, it does seem clear 
that three-dimensional effects such as the formation of streamwise ligaments and 
the pulling back of the sheet at its edges under the action of surface tension are 
important. 

A large body of analytic work on the linear stability of plane sheets exists (Squire, 
1953; Hagerty and Shea, 1955). These papers found the sinuous mode (antisym- 
metric mode) to be most unstable but concluded that disintegration of the sheet 
could not be predicted. This follows since the linear analysis does not predict thin- 
ning of the sheet. The nonlinear planar problem has been studied using asymptotic 
methods (Clark and Dombrowski, 1972), and more recently a reduced dimension 
technique has been used by, Mehring and Sirigano (1999). These analyses predict a 
thinning of the sheet that may lead to the sheet breakup and disintegration. Anal- 
ysis of more general geometries has been limited although Mehring and Sirigano 
(2000) have recently started work on analysis of a conical sheet with their reduced- 
dimension technique. In the future, analytic work with directed fluid sheets using 
techniques similar to those used by Bogy (1979) for liquid jets may be of some use. 

Modeling the expected behavior of the conical sheet will require a fully three- 
dimensional code. Numerical simulations of the breakup of three-dimensional liq- 
uid sheets have been limited. These calculations are difficult due to problems with 
accurately tracking fluid interfaces and gross topological changes. Several tech- 
niques have evolved to address the difficulties associated with having a surfaces 
of discontinuity embedded in a three-dimensional domain. Such methods include: 
Level Set methods, Sethian (1999), Volume of Fluid methods Zaleski (1999) and 
Front Tracking schemes, Unverdi & Tryggvason (1992). Finite element techniques 
have also been used for interface problems and offer attractive accuracy but aie 
computationally very expensive, and, so far, topological changes have not been 
implemented. 

For our work we have chosen to develop a front tracking code based on the method 
used by Tryggvason’s group. This method tracks material particles on the interface 
from which material property fields such as density and viscosity are constructed 
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Fig. 1. A fixed two dimensional grid with an embedded surface constructed from 
material points separating two regions of fluid with disparate material properties. 

at each time step. A field of surface forces is also calculated based on the front’s 
location and geometry. 

By keeping track of material points on the interface, this scheme also offers a 
means to control topological changes and introduce physics that control such occur- 
rences. These changes are typically due to processes at the molecular level. Adding 
this capability to front tracking schemes is currently an active area of research. 

We are also interested in developing this code to understand the instability of 
a thin liquid film lying between a solid surface and an incompressible gas flowing 
over it. This has applications in atomizers and hybrid rocket engines. For this 
situation we intend to use a planar two-dimensional simulation to understand the 
instability and predict droplet formation. The results given below are from the 
interface tracking scheme to be used with this code. 

2. Accomplishments 

In this section we describe the two-dimensional code we are writing using the 
interface tracking technique developed by Tryggvason’s group (Uneverdi and Tryg- 
gvason, 1992). Although this section describes a two-dimensional interface tracking 
method, it illustrates the basic methodology used for the three-dimensional case. 

A unique feature of interface tracking schemes is the ability to follow the motion 
of material particles embedded in the interface as shown in Fig. 1. This introduces 
several extra complexities into the code needed to keep track of the location and 
properties of the surface. However, language features such as pointers and structures 
found in modern programming languages allow us to easily construct data structures 
that greatly simplify this task. 

2.1 Front tracking algorithm 

The front tracking scheme uses a fixed grid through which a material interface 
is convected. The material interface is represented by a discrete series of points 
representing the location of the front. Material properties are allowed to jump 
across the interface. The material particles are convected at each time step, and 
new material property fields and surface force fields are calculated based on the 
position of points on the front. 
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Equations of motion 

A form of the Navier- Stokes equations that includes effects of stresses intrinsic to 
a surface of discontinuity embedded in the domain is written as 


<3 u 

— h div(pu 0 u) = — gradp + pi + div[p(gradu + gradu T )] 
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Here p represents density, p viscosity, p pressure, u fluid velocity, and f are body 
forces. 7 is surface tension, H is mean curvature, and r is the position vector of an 
arbitrary point in the domain while r 3 is a point on the surface. The vector n is the 
unit normal to the surface of the interface, and is the multi-dimensional delta 
function. This equation is valid over the entire domain (Unverdi & Tryggvason 
1992). 

Material properties of material particles are assumed to remain constant in time 
so that p = 0 and p — 0. Now the mass conservation equation 


^ + div(pu) = 0 , 


simply becomes 

div(u) = 0 . ( 2 ) 

Equation 1 is discretized onto the fixed grid using standard schemes. An impor- 
tant difference with standard incompressible DNS codes is that the effects of surface 
forces are included on the fixed grid and material properties are spatially variable 
although the flow is incompressible. The density variation leads to a non-separable 
PDE for the pressure equation that does not allow the use of fast Poisson-type 
solvers. As a result, the pressure equation is usually solved using multi-grid solvers 
such as MUDPAK, which substantially increases computational time. 

Description of the front 

The front is represented by a series of FORTRAN 90, TYPE structures connected 
by pointers to form a doubly-linked list. A linked list is a natural metaphor for the 
material front; each node of the structure contains information pertinent to that 
point on the interface. Variables such as position, arc-length, and varying intrinsic 
surface properties can easily be stored in each node structure. 

Each node of the linked list represents an endpoint of a surface element. This 
data structure allows us to easily insert or remove surface elements by inserting 
nodes as the front is stretched or compressed. The insertion-extraction scheme is 
adjusted to give about four elements per grid-cell. The position of new nodes is 
determined by simple polynomial interpolation. 

Construction of material property fields 

Fluid properties in the bulk regions are constructed at each time step based on 
the updated location of the interface. This includes constructing a field of surface 
forces that are transferred from the front to the fixed grid. 
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The scheme for constructing the density fiehl follows. First, we note that the 
density gradient can be expressed as 

gradp = j Apn6(r — r 3 )ds 

where A p is the jump in density across the interface. The subscript s denotes values 
on the surface, n is the unit normal as used in Eq. 1 and is oriented such that 

gradp 

h = ligradpjf' 

A value for gradp on the fixed grid is calculated by using a suitable weighting 
function, uq ; *, to transfer information from the front to the fixed grid. Wijk is the 
product of n one-dimensional functions: 


w, } k(r 3 ) ~ d(x 3 - x)d(y 3 - y)d{z 3 - z). 


(3) 


Here x, y, z are the coordinates of the fixed grid and x 3 , y s , z 3 are coordinates of a 
point on the surface. r 3 is the position of a point on the surface so that Wijk( r 3 ) is 
the weighting on a fixed grid point due to a surface element at r 3 . In our work so 
far we have used the function first used by Peskin (1977) and now commonly used: 


d{r) 


( l + cos(nr/2h ) 
2 h 

o, 


M < 2h , 

jjrjj > 2 h, 


(4) 


where r = ||r — r s ||. 

After visiting each node on the interface, we are left with a smoothed approx- 
imation to the derivative across the interface. Examples of the smoothed density 
gradient are shown in Fig. 2 and Fig. 3. The density ratio is 10:1 with the fluid 
surrounding the drop assumed to have a unit density. 

Taking the numerical divergence of the density gradient yields a Poisson equation 
for the density field: 

div(gradp) = div^grad p i} (5) 

The term on the right denotes the numerical divergence of the grid density gradient 
field. This equation is readily solved for p with appropriate boundary conditions 
using fast-Poisson solvers such as the routines found in FISHPACK. The smoothed 
reconstructed field is shown in Fig. 4 . The viscosity field is reconstructed in the 
same way. 

Surface forces 

Numerical simulations of interfacial problems with surface tension introduce sev- 
eral difficulties. In the scheme discussed here surface forces are transferred to the 
front in a similar manner as discussed above for the density field. 


i 
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Fig. 5. Surface forces transferred to the fixed 64 x 64 grid from an initially 
circular drop. 



Fig. 6. Surface forces transferred to a fixed grid after the drop has been deformed 
in a stagnation point flow. Several surface points have been added in areas where 
the surface has been stretched and several deleted in the regions of higher curvature 
where the surface has been compressed. 
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For two-dimensional in the plane, mean curvature simply becomes the curvature 
of the plane curve marking the interface. The force transferred to the bulk fluid 
due to a surface element of length ds is then given by 

j\ lAlU = ,I^JM ds= , {t2 . h) . (6, 

Here, k is the curvature and n is the normal to the plane curve in the sense of 
the Frenet formulae. That is, it points toward the center of curvature of the plane 
curve and not necessarily in the direction of n, which acquires its direction from a 
gradient in the density field. 

As with the density field, we visit each surface element in the linked list that 
represents the front. At each node, the unit tangents tj and t 2 at each end of the 
line element are calculated. In two dimensions, the integrated effect of the surface 
tension forces over the element length is simply the vector difference in the unit 
tangent vectors of node points at the end of the surface element as shown in Eq. 6. 

After visiting all the points on the front we have a field of surface forces transferred 
to the fixed grid as shown in Fig. 5 and Fig. 6. Figure 5 illustrates the surface force 
field created by an initially circular “drop”. Figure 6 shows the surface force field 
after the drop has been strained into an ellipse. The greater values of surface forces 
transferred to the fixed grid are apparent in the deformed case due to the higher 
curvatures at the right and left ends of the drop. 

3. Future plans 

In the next year we intend to complete a two-dimensional planar and axisymmet- 
ric front-tracking code. This code will be used to study the instability and drop 
formation of a thin sheet lying on a solid surface. This work will expand on the 
present effort, leading to more sophisticated models that incorporate phase change 
as well as topological changes. 

A substantial challenge for the forthcoming year will be the implementation of 
a fully three-dimensional code, which is ultimately what we require to study the 
conical jet. As with the two-dimensional code, the greatest obstacle to overcome 
is accurate representation of the front and construction of fluid property fields. 
In three-dimensions the primary difference is that the front is represented as a 
collection of two-dimensional surface elements. The size and distribution of these 
elements must be adjusted with time as surface elements become too small or large. 

A major difficulty in three-dimensional codes is the accurate calculation of surface 
curvature. At present most front- tracking schemes use fitted ellipsoids to estimate 
the local curvature. We hope to construct a scheme that uses the coordinates of 
the centers of the surface elements. 

Also, alternative techniques have been developed by Popinet and Zaleski (1999) 
to transfer surface forces onto the fixed grid in two-dimensions. These authors 
have reported that this technique reduces spurious currents found in other schemes 
due to smoothing surface, stresses onto the fixed grid. Their techniques should, in 
principle, be adaptable to three-dimensional cases. 
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In the front tracking scheme, the underlying differencing grid is a uniform station- 
ary mesh. In general, to obtain a manageable computation, the grid is a compromise 
between resolving the small length scales at the interface and computational power. 
Adapting the grid to a finer resolution near the interface where relevant length 
scales are small compared to the bulk flow is one solution. This technique has been 
implemented in two dimensions by Agresar et a l. (1998) to study biological cell 
mechanics. Its implementation in three-dimensions will be investigated. 
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